library(tidyverse)
library(gtable)
library(gtsummary)
library(gridExtra)
library(knitr)
library(car)
library(boot)
data <- readRDS("Dataset.RData")

melanoma <- readRDS("Dataset.RData") %>%
  mutate(age_c = mean(age)-age,
         trans_class_bin = ifelse(trans_class == "immune",1,0),
         prop_infiltrated = infiltration_count/tile_count,
         stage_bin = ifelse(stage == 3, "3",
                            ifelse(stage == 4, NA, "1-2"))) %>%
  filter(!is.na(stage_bin))

data <- melanoma

N = length(melanoma$id)
data$y = 0
data$y[data$survival == 'Long'] = 1
# immune selected genes MC1R, AKT1, NAMPT, TNFRSF19, PAK2, SOS1, B2M 
l1 <- glm(data$y ~ data$MC1R + as.factor(data$stage), family = 'binomial')
l1$coefficients[2]
##  data$MC1R 
## -0.9632995
l2 <- glm(data$y ~ data$AKT1 + as.factor(data$stage), family = 'binomial')
l2$coefficients[2]
## data$AKT1 
## -0.917543
l3 <- glm(data$y ~ data$NAMPT + as.factor(data$stage), family = 'binomial')
l3$coefficients[2]
## data$NAMPT 
##  0.1939064
l4 <- glm(data$y ~ data$TNFRSF19 + as.factor(data$stage), family = 'binomial')
l4$coefficients[2]
## data$TNFRSF19 
##       0.23159
l5 <- glm(data$y ~ data$PAK2 + as.factor(data$stage), family = 'binomial')
l5$coefficients[2]
## data$PAK2 
## 0.5984475
l6 <- glm(data$y ~ data$SOS1 + as.factor(data$stage), family = 'binomial')
l6$coefficients[2]
## data$SOS1 
## 0.2604371
l7 <- glm(data$y ~ data$B2M + as.factor(data$stage), family = 'binomial')
l7$coefficients[2]
##  data$B2M 
## 0.3445805
betas <- unname(c(l1$coefficients[2], l2$coefficients[2], l3$coefficients[2], l4$coefficients[2], l5$coefficients[2], l6$coefficients[2], l7$coefficients[2]))
dd = data[c('MC1R', 'AKT1', 'NAMPT', 'TNFRSF19', 'PAK2', 'SOS1', 'B2M')] 
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(ggcorrplot)
library(heatmaply)
## Loading required package: viridis
## Loading required package: viridisLite
## 
## ======================
## Welcome to heatmaply version 1.3.0
## 
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
## 
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags: 
##   https://stackoverflow.com/questions/tagged/heatmaply
## ======================
heatmaply_cor(cor(dd))
cor(dd)
##                MC1R       AKT1      NAMPT   TNFRSF19       PAK2       SOS1
## MC1R      1.0000000  0.6168210 -0.1632614 -0.1315617 -0.3218287 -0.3567775
## AKT1      0.6168210  1.0000000 -0.2856350 -0.2455591 -0.4477137 -0.3738679
## NAMPT    -0.1632614 -0.2856350  1.0000000  0.2330692  0.1857005  0.1142098
## TNFRSF19 -0.1315617 -0.2455591  0.2330692  1.0000000  0.2732336  0.0416091
## PAK2     -0.3218287 -0.4477137  0.1857005  0.2732336  1.0000000  0.2164186
## SOS1     -0.3567775 -0.3738679  0.1142098  0.0416091  0.2164186  1.0000000
## B2M      -0.2274429 -0.1892508  0.2003606 -0.0800356  0.1264879  0.1687312
##                 B2M
## MC1R     -0.2274429
## AKT1     -0.1892508
## NAMPT     0.2003606
## TNFRSF19 -0.0800356
## PAK2      0.1264879
## SOS1      0.1687312
## B2M       1.0000000
dd1 = as.data.frame(dd)
dd1$y = data$y
dd1$age_c = data$age - mean(data$age)
dd1$sex = 0
dd1$sex[data$gender == 'female'] = 1
dd1$pathStage = as.factor(data$stage_bin)
model = glm(y ~ MC1R + AKT1 + NAMPT + TNFRSF19 + PAK2 + SOS1 + B2M + pathStage, data = dd1, family = 'binomial')
data.frame(vif(model)[-8], betas) %>%
  kable(col.names = c("VIF","$\\hat{\\beta}$"))
VIF \(\hat{\beta}\)
MC1R 1.182171 -0.9632995
AKT1 1.390206 -0.9175430
NAMPT 1.127706 0.1939064
TNFRSF19 1.189646 0.2315900
PAK2 1.183120 0.5984475
SOS1 1.183708 0.2604371
B2M 1.098849 0.3445805
# selected top genes: MC1R, AKT1, PAK2, SOS1, B2M
# model 1: basline + selected top genes

model1 <- glm(y ~ MC1R + AKT1 + PAK2 + sex + age_c + pathStage, data = dd1, family = 'binomial')
summary(model1)
## 
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + sex + age_c + pathStage, 
##     family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2138  -0.9441   0.4158   0.8805   1.9620  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.68720    0.24066   2.856 0.004297 ** 
## MC1R        -0.53926    0.28790  -1.873 0.061055 .  
## AKT1        -0.56258    0.23008  -2.445 0.014481 *  
## PAK2         0.40056    0.17698   2.263 0.023617 *  
## sex         -0.22564    0.30450  -0.741 0.458669    
## age_c       -0.02940    0.01021  -2.880 0.003977 ** 
## pathStage3  -1.18109    0.31076  -3.801 0.000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 330.38  on 238  degrees of freedom
## Residual deviance: 265.65  on 232  degrees of freedom
## AIC: 279.65
## 
## Number of Fisher Scoring iterations: 5
model1.diag <- glm.diag(model1)
data.frame("cooks" = model1.diag$cook, "leverage" = model1.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/N)) %>%
  ggplot(aes(x = seq(0.5,N,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

N=239
p_res1 <- residuals(model1, type = "pearson")
d_res1 <- residuals(model1, type = "deviance")
data.frame("pearson" = p_res1, "deviance" = d_res1) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model1.diag$rp, "deviance" = model1.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res1, "deviance" = d_res1,"leverage" = model1.diag$h, "cooks" = model1.diag$cook,"p_standard" = model1.diag$rp, "d_standard" = model1.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 2: baseline + selected top genes + immune response genes
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$DCK = data$DCK
dd1$C5 = data$C5
dd1$CCL28 = data$CCL28
dd1$RFX5 = data$RFX5
dd1$EIF2AK2 = data$EIF2AK2
model2 <- glm(y ~ MC1R + AKT1 + PAK2 + sex + age_c + DDX58 + KLRD1 + TLR4 + pathStage, data = dd1, family = 'binomial')
summary(model2)
## 
## Call:
## glm(formula = y ~ MC1R + AKT1 + PAK2 + sex + age_c + DDX58 + 
##     KLRD1 + TLR4 + pathStage, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1888  -0.9154   0.3934   0.8375   2.1739  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.82951    0.25211   3.290  0.00100 ** 
## MC1R        -0.27151    0.29388  -0.924  0.35554    
## AKT1        -0.55906    0.23691  -2.360  0.01829 *  
## PAK2         0.32341    0.18538   1.745  0.08105 .  
## sex         -0.23917    0.31225  -0.766  0.44370    
## age_c       -0.03071    0.01056  -2.908  0.00364 ** 
## DDX58        0.04621    0.16474   0.281  0.77908    
## KLRD1        0.26288    0.18681   1.407  0.15938    
## TLR4         0.38050    0.21225   1.793  0.07302 .  
## pathStage3  -1.38949    0.33121  -4.195 2.73e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 330.38  on 238  degrees of freedom
## Residual deviance: 256.07  on 229  degrees of freedom
## AIC: 276.07
## 
## Number of Fisher Scoring iterations: 5
model2.diag <- glm.diag(model2)
data.frame("cooks" = model2.diag$cook, "leverage" = model2.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
  ggplot(aes(x = seq(0.5,N,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

p_res2 <- residuals(model2, type = "pearson")
d_res2 <- residuals(model2, type = "deviance")
data.frame("pearson" = p_res2, "deviance" = d_res2) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model2.diag$rp, "deviance" = model2.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res2, "deviance" = d_res2,"leverage" = model2.diag$h, "cooks" = model2.diag$cook,"p_standard" = model2.diag$rp, "d_standard" = model2.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
# model 3: baseline + selected top genes (2!!!) + immune response genes + age:immune response genes
model3 <- glm(y ~ pathStage*MC1R + pathStage*AKT1 + pathStage*PAK2 + sex + age_c + pathStage*DDX58 + pathStage*KLRD1 + pathStage*TLR4, data = dd1, family = 'binomial')

model4<- glm(y ~ pathStage*MC1R + pathStage*AKT1 + pathStage*PAK2 + sex + age_c, data = dd1, family = 'binomial')

summary(model3)
## 
## Call:
## glm(formula = y ~ pathStage * MC1R + pathStage * AKT1 + pathStage * 
##     PAK2 + sex + age_c + pathStage * DDX58 + pathStage * KLRD1 + 
##     pathStage * TLR4, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9111  -0.8781   0.1851   0.8371   2.0644  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       1.06409    0.30921   3.441 0.000579 ***
## pathStage3       -1.57509    0.39032  -4.035 5.45e-05 ***
## MC1R             -0.18903    0.37501  -0.504 0.614220    
## AKT1             -0.54072    0.32124  -1.683 0.092331 .  
## PAK2              0.51523    0.28225   1.825 0.067931 .  
## sex              -0.12842    0.32561  -0.394 0.693273    
## age_c            -0.03258    0.01102  -2.958 0.003100 ** 
## DDX58            -0.27738    0.20877  -1.329 0.183966    
## KLRD1             0.42886    0.44573   0.962 0.335973    
## TLR4              0.95109    0.48074   1.978 0.047884 *  
## pathStage3:MC1R  -0.10106    0.60318  -0.168 0.866944    
## pathStage3:AKT1  -0.05622    0.48646  -0.116 0.907988    
## pathStage3:PAK2  -0.51131    0.39054  -1.309 0.190461    
## pathStage3:DDX58  1.08710    0.45414   2.394 0.016678 *  
## pathStage3:KLRD1 -0.23974    0.49227  -0.487 0.626260    
## pathStage3:TLR4  -0.98786    0.55348  -1.785 0.074290 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 330.38  on 238  degrees of freedom
## Residual deviance: 244.71  on 223  degrees of freedom
## AIC: 276.71
## 
## Number of Fisher Scoring iterations: 6
summary(model4)
## 
## Call:
## glm(formula = y ~ pathStage * MC1R + pathStage * AKT1 + pathStage * 
##     PAK2 + sex + age_c, family = "binomial", data = dd1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3949  -0.9524   0.3473   0.8835   1.9123  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.71413    0.24623   2.900 0.003728 ** 
## pathStage3      -1.21018    0.33113  -3.655 0.000258 ***
## MC1R            -0.49375    0.36101  -1.368 0.171416    
## AKT1            -0.51748    0.30760  -1.682 0.092506 .  
## PAK2             0.62292    0.26121   2.385 0.017089 *  
## sex             -0.20928    0.30744  -0.681 0.496045    
## age_c           -0.02825    0.01030  -2.742 0.006099 ** 
## pathStage3:MC1R -0.02713    0.59157  -0.046 0.963426    
## pathStage3:AKT1 -0.11487    0.46470  -0.247 0.804758    
## pathStage3:PAK2 -0.43708    0.36093  -1.211 0.225905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 330.38  on 238  degrees of freedom
## Residual deviance: 264.15  on 229  degrees of freedom
## AIC: 284.15
## 
## Number of Fisher Scoring iterations: 6
anova(model3, model4, test = "LRT")
anova(model1, model2, test = "LRT")
model3.diag <- glm.diag(model3)
data.frame("cooks" = model3.diag$cook, "leverage" = model3.diag$h) %>%
  pivot_longer(everything(), names_to = "measure") %>%
  mutate(cutoff = ifelse(measure=="cooks",1,26/240)) %>%
  ggplot(aes(x = seq(0.5,N,by=0.5))) +
  geom_line(aes(y = value, col = measure)) +
  geom_line(aes(y = cutoff), lty = 2, alpha = 0.5) +
  theme_classic() +
  facet_wrap(~measure, scales = "free_y", nrow = 2) +
  scale_color_manual(values = c("orange","cornflowerblue")) +
  labs(x = "observation")

p_res3 <- residuals(model3, type = "pearson")
d_res3 <- residuals(model3, type = "deviance")
data.frame("pearson" = p_res3, "deviance" = d_res3) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen","purple")) +
  labs(x = "observation")

data.frame("pearson" = model3.diag$rp, "deviance" = model3.diag$rd) %>%
  pivot_longer(everything(), names_to = "residual") %>%
  ggplot(aes(x = seq(0.5,N,by=0.5), y = value, col = residual)) +
  geom_point() +
  facet_wrap(~residual, nrow = 2) +
  theme_classic() +
  scale_color_manual(values = c("darkgoldenrod2","orchid3")) +
  labs(x = "observation", y = "standardized values")

#kable(data.frame("pearson" = p_res3, "deviance" = d_res3,"leverage" = model3.diag$h, "cooks" = model3.diag$cook,"p_standard" = model3.diag$rp, "d_standard" = model3.diag$rd),col.names = c("$\\hat{r}^{P}$","$\\hat{r}^{D}$","leverage","cooks","$\\hat{r}^{PS}$","$\\hat{r}^{DS}$"))
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following objects are masked from 'package:boot':
## 
##     logit, simplex
## The following object is masked from 'package:car':
## 
##     logit
## The following object is masked from 'package:tidyr':
## 
##     fill
# model 5: proportional odds model
dd1$CD86 = data$CD86
dd1$DDX58 = data$DDX58
dd1$KLRD1 = data$KLRD1
dd1$TLR4 = data$TLR4
dd1$pathStage = as.factor(data$stage)
model5 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, data = dd1, family = cumulative(parallel=T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model5)
## 
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + 
##     sex + age_c + CD86 + DDX58 + KLRD1 + TLR4, family = cumulative(parallel = T), 
##     data = dd1)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.196474   0.156455  -7.647 2.05e-14 ***
## (Intercept):2  0.307711   0.139717   2.202   0.0276 *  
## MC1R          -0.315414   0.144729  -2.179   0.0293 *  
## AKT1           0.285869   0.125643   2.275   0.0229 *  
## PAK2           0.249782   0.105872   2.359   0.0183 *  
## SOS1           0.065107   0.099621   0.654   0.5134    
## B2M            0.095804   0.131593   0.728   0.4666    
## sex           -0.204166   0.186470  -1.095   0.2736    
## age_c         -0.012211   0.005872  -2.079   0.0376 *  
## CD86          -0.323995   0.164611  -1.968   0.0490 *  
## DDX58          0.234769   0.100119   2.345   0.0190 *  
## KLRD1          0.031465   0.138238   0.228   0.8199    
## TLR4          -0.058970   0.117485  -0.502   0.6157    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2])
## 
## Residual deviance: 489.8513 on 465 degrees of freedom
## 
## Log-likelihood: -244.9256 on 465 degrees of freedom
## 
## Number of Fisher scoring iterations: 14 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##      MC1R      AKT1      PAK2      SOS1       B2M       sex     age_c      CD86 
## 0.7294869 1.3309179 1.2837450 1.0672736 1.1005435 0.8153270 0.9878635 0.7232542 
##     DDX58     KLRD1      TLR4 
## 1.2646165 1.0319651 0.9427350
model6 <- vglm(pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + sex + age_c + CD86 + DDX58 + KLRD1 + TLR4
                + age_c:CD86 + age_c:DDX58 + age_c:KLRD1 + age_c:TLR4 , data = dd1, family = cumulative(parallel = T))
## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()
summary(model6)
## 
## Call:
## vglm(formula = pathStage ~ MC1R + AKT1 + PAK2 + SOS1 + B2M + 
##     sex + age_c + CD86 + DDX58 + KLRD1 + TLR4 + age_c:CD86 + 
##     age_c:DDX58 + age_c:KLRD1 + age_c:TLR4, family = cumulative(parallel = T), 
##     data = dd1)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.227502   0.164555  -7.460 8.68e-14 ***
## (Intercept):2  0.351872   0.146671   2.399 0.016437 *  
## MC1R          -0.378942   0.146742  -2.582 0.009812 ** 
## AKT1           0.277384   0.130332   2.128 0.033313 *  
## PAK2           0.249584   0.109425   2.281 0.022556 *  
## SOS1           0.071908   0.103761   0.693 0.488302    
## B2M            0.137753   0.141915   0.971 0.331713    
## sex           -0.251573   0.195706  -1.285 0.198630    
## age_c         -0.019616   0.006452  -3.040 0.002365 ** 
## CD86          -0.233614   0.171843  -1.359 0.173999    
## DDX58          0.329202   0.117563   2.800 0.005107 ** 
## KLRD1         -0.091329   0.162974  -0.560 0.575212    
## TLR4          -0.113594   0.126943  -0.895 0.370871    
## age_c:CD86     0.021738   0.013329   1.631 0.102925    
## age_c:DDX58   -0.035834   0.010462  -3.425 0.000614 ***
## age_c:KLRD1   -0.028273   0.010840  -2.608 0.009104 ** 
## age_c:TLR4    -0.011448   0.011069  -1.034 0.301021    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2])
## 
## Residual deviance: 474.2647 on 461 degrees of freedom
## 
## Log-likelihood: -237.1323 on 461 degrees of freedom
## 
## Number of Fisher scoring iterations: 15 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##        MC1R        AKT1        PAK2        SOS1         B2M         sex 
##   0.6845854   1.3196733   1.2834914   1.0745560   1.1476919   0.7775766 
##       age_c        CD86       DDX58       KLRD1        TLR4  age_c:CD86 
##   0.9805748   0.7916670   1.3898589   0.9127171   0.8926203   1.0219761 
## age_c:DDX58 age_c:KLRD1  age_c:TLR4 
##   0.9648000   0.9721226   0.9886169